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Abstract 

Highly  accurate  finite  difference  schemes  are  developed  for 
Laplace's  equation  with  the  Dirichlet  boundary  condition  on  general 
bounded  regions  in  R  .   A  second  order  accurate  scheme  is  combined 
with  a  deferred  correction  or  Richardson  extrapolation  method  to 
increase  the  accuracy.   The  Dirichlet  condition  is  approximated  by 
a  method  suggested  by  Heinz-Otto  Kreiss.   A  convergence  proof  of 
his,  previously  not  published,  is  given  which  shows  that,  for  the 
interval  size  h,  one  of  the  methods  has  an  accuracy  of  at  least 
0{yP      )    in  Lo.   The  linear  systems  of  algebraic  equations  are 
solved  by  a  capacitance  matrix  method.   The  results  of  our  numeri- 
cal experiments  show  that  highly  accurate  solutions  are  obtained 
with  only  a  slight  additional  use  of  computer  time  when  compared 
to  the  results  obtained  by  second  order  accurate  methods. 


AMS  (MOS)  subject  classifications.   Primary  65BO5,  65NI5; 
secondary  65FO5,  65N20. 


§1.  Introduction 

It  is  the  purpose  of  this  paper  to  develop  some  highly 
accurate  finite  difference  methods  for  the  Dirichlet  problem  for  a 
general  bounded  region  O  in  R  .   The  most  accurate  of  these  has 
an  L^  error  of  order  at  most  h^'^,  see  §4.   Our  basic  schemes  use 
the  standard  (2n+l)-point  formula  for  the  interior  mesh  points  and 
are  therefore  only  second  order  accurate.   The  increased  accuracy 
is  achieved  by  two  steps  of  a  deferred  correction  or  Richardson 
extrapolation  procedure.   We  also  discuss  the  computer  implementa- 
tion of  these  methods  in  some  detail. 

The  use  of  deferred  correction  and  Richardson  extrapolation 
methods  is  justified  by  finding  asymptotic  expansions  of  the  error. 
Wasow  [20]  has  shown  that  no  useful  expansions  of  this  kind  exist 
if  the  boundary  condition  is  approximated  to  a  low  order  of 
accuracy.   An  obvious  remedy  for  this  problem,  already  mentioned  by 
Wasow,  is  to  use  higher  order  interpolation  or  extrapolation  formu- 
las at  any  irregular  mesh  point,  i.e.  a  mesh  point  in  the  open  set 
O  which  fails  to  have  all  its  next  neighbors  in  the  closure  of  O. 
Volkov  [19]  proposed  the  use  of  high  order  one-dimensional  Lagrange 
polynomials  for  this  purpose.   Because  of  the  change  of  sign  of  the 
interpolation  coefficients  the  matrix  representing  the  difference 
scheme  will  then,  in  general,  not  be  of  positive  type.   The  stan- 
dard convergence  proof  based  on  a  discrete  maximum  principle,  see 
Forsythe  and  Wasow  [7],  will  therefore  generally  not  apply.   But  by 
allowing  the  use  of  values  of  the  mesh  functions  many  mesh  lengths 
away  from  the  boundary,  Volkov  succeeded  in  designing  schemes  with 


diagonally  dominant  matrices.   His  schemes  may  however  lead  to  an 
unacceptably  small  mesh  size  even  for  very  simple  geometries. 

Numerical  experiments,  see  Pereyra  [13]    and  the  last 
section  of  this  paper,  clearly  demonstrate  the  need  for  higher 
order  accuracy  at  the  irregular  mesh  points  if  improved  solutions 
through  Richardson  extrapolation  or  deferred  correction  methods 
are  required.   In  his  1966  paper,  Pereyra  also  reported  on 
successful  numerical  experiments  with  methods  based  on  Lagrange 
interpolation  in  one  variable  and  employing  only  mesh  points  close 
to  the  boundary.   At  that  time  no  convergence  proof  was  known  for 
such  methods. 

In  June  of  1968,  Kreiss  announced  an  interesting  result  on 
the  convergence  of  methods  of  this  type.   His  result  was  never 
published.   His  schemes  are  constructed  as  sums  of  difference 
approximations  of  one-dimensional  problems.   At  the  interior  mesh 
points  each  of  these  problems  is  discretized  by  a  three-point 
formula  while  at  the  irregular  mesh  points  this  basic  formula  is 
combined  with  high  order  Lagrange  extrapolation.   For  a  detailed 
description  see  §2.   Kreiss  found  a  method  of  proof  which  provides 
an  alternative  to  the  classical  technique  previously  mentioned. 
His  method  depends  heavily  on  the  special  structure  just  described. 

We  learned  about  his  results  from  several  conversations  and 
his  unpublished  notes  which  were  kindly  made  available  to  us.   Our 
interest  in  these  methods  was  recently  renewed  when  we  realized 
that  the  capacitance  matrix,  or  imbedding,  method  developed  by 
Proskurowski  and  Widlund  [18]  could  be  adapted  for  the  difference 
schemes  considered  by  Kreiss. 


In  this  paper,  we  describe  Kreiss'  schemes,  give  detailed 
proofs  of  convergence  and  existence  of  error  expansions  and 
discuss  their  implementation.   We  have  exclusively  used  a  deferred 
correction  method  in  our  numerical  experiments  rather  than  Richardson 
extrapolationo   Our  reason  is  that  the  deferred  correction  method, 
especially  for  problems  in  several  dimensions,  has  often  proved 
less  costly,  see  Pereyra  [13]  and  also  §  5  of  this  paper.   One 
advantage  is  that,  in  contrast  to  Richardson  extrapolation, 
deferred  correction  methods  require  only  one  mesh  size.   The 
capacitance  matrix  method  allows  us  to  solve  the  same  system  of 
linear  equations  repeatedly  at  an  expense  which  decreases  con- 
siderably once  the  first  problem  has  been  solved. 

Our  combination  of  a  deferred  correction  and  an  imbedding 
method  is  quite  convenient  from  a  programming  point  of  view.   We 
have  also  developed  a  new,  practical  way  of  calculating  the 
required  difference  approximations  to  the  terms  of  the  expansion 
of  the  truncation  error.   This  method  resolves  a  long-standing 
problem  in  the  theoretical  justification  for  the  use  of  more  than 
one  deferred  correction  step  for  boundary  value  problems  of  this 
type.   The  imbedding  of  the  region  in  a  rectangle  allows  us  to  use 
certain  programs  previously  developed  to  perform  deferred  correc- 
tions for  problems  on  rectangular  regions. 

In  the  last  section, we  report  on  numerical  experiments 
carried  out  on  a  CDC  760O  computer  at  the  Lawrence  Berkeley 
Laboratory.   They  show  that  very  high  accuracy  is  obtained  for 
problems  with  sufficiently  smooth  solutions.   For  problems  which 
fail  to  have  sufficiently  many  bounded  derivatives  the  corrections 


do  not  spoil  the  accuracy  of  the  solution.   We  believe  that  our 
method  can  "be  developed  further  into  highly  efficient  and  reliable 
numerical  software.   We  note  that  fast  Laplace  solvers  are  used 
increasingly  to  enhance  the  convergence  when  solving  more  general 
problems,  see  for  example,  Bartels  and  Daniel  [1],  Concus  and  Golub 
[4,5],  Jameson  [9],  O'Leary  [10],  Martin  [11,12]  and  Widlund  [21]. 
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§2o  Kreiss'  Method  for  Poisson's  Equation 

We  will  consider  a  family  of  finite  difference  schemes  for 
the  Dirichlet  problem  for  Poisson's  equation. 


_Au  =  -  J~  ih/hx^)    u  =  f(x)  ,    X  e  O  , 
(2.1) 

u(x)  =  g(x)  ,    X  e  ^  , 

where  the  region  O  is  an  open,  bounded  subset  of  the  n-dimensional, 
real,  Euclidean  space  R  with  the  boundary  ^.   We  will  make  no 
detailed  assumptions  on  the  smoothness  of  ^O  and  the  data  f  and  g 
but  assume  only  that  they  are  sufficiently  smooth.   As  is  well 
known,  the  problem  (2.1)  then  has  a  unique,  sufficiently  smooth 
solution, 

A  uniform  mesh  R,  is  introduced  by 

r!!?  =  fx  e  r'^I  X.  -  n.h,  n.  =  0,±1,±2,...]  , 

where  h  >  0  is  the  mesh  size.   The  position  of  the  origin  of  our 
mesh  is,  of  course,  arbitrary.   We  could  also  have  chosen  different 
uniform  mesh  sizes  in  the  different  coordinate  directions  without 
affecting  the  theory  or  practice  of  the  methods  except  in  some  very 
minor  ways. 

The  set  of  mesh  points  of  interest  to  us  is 


n 

=  r  1  I  I  K 


o^  =  o  nRj^ 


There  are  no  equations  for  points  on  hO,      The  difference  equations 
are  constructed  as  a  sum  of  approximations  of  one-dimensional  prob- 


lems  corresponding  to  the  operators  -(d/Sx.)  ,  j=l,  ...,n.  They  are 
specified,  by  defining  a  linear  equation  for  each  x  e  Oy^.   Let  the 
vector  e.  "be  the  unit  vector  in  the  direction  of  the  positive  i-th 


closest  neighbors  x±he.,  i  =  l,...,n,  belong  to  O^.   For  a  regular 
mesh  point,  we  simply  use  the  standard  centered  difference 
approximation  of  each  of  the  second  derivatives.   This  results  in 
the  equation 

2nu^(x)  -J—   (u^(x+he^)+  u^(x-  he^))  =h^f(x)  . 

This  formula  is  combined  with  polynomial  extrapolation  of  a 

fixed  degree  k  for  the  remaining,  irregular,  mesh  points  of  O,  o 

Let  us  thus  suppose  that  x  e  O,  but  that  x+  he.  /  O,  and  that 

X  -  he  . ,  . .  o ,  X  -  (k-l)he.  e  O,  .   This  last  condition  can  8.1ways  be 

satisfied  for  a  smooth  hO   if  h  is  chosen  small  enough.   Denote  by 

X?  the  intersection  of  the  boundary  hO   and  the  segment  between  x 

and  x+he.    and  by     s.h    the    distance   between  x*   and  x  +  he .  .      Thus 
1  "^  11 

0  ;<   s   <   1.      A  provisional  value   of   u    (x  +  he^)    is   now   defined  by   the 
Lagrange    interpolation   formula, 

(2.2)  ^  a   u^(x-(j-l)he.)    =   u(x*)=g(x*)    . 

The  coefficients  a  .  depend  cnly  on  s  and  are  given  by  the  formula 

k 
a.  =   TT  (s-^)/(j-i)  . 

^    £=0 

The  value  of  u^  at  the  point  x+  he.  is  now  eliminated  by  combining 


(2.2)  with  the  standard  three  point  formula  for  the  point  x.   The 
resulting  matrix,  which  corresponds  to  the  approximation  of 

o 

-{h/hx.)      along  a  mesh  line  parallel  to  e.,thus  typically  has  the 
form 

(2.3) 


(2+a^/a^),  (-l+a2/a^),a^/a^,  ...,a^/a^ 


-1 


Q^/a^,  ...,a^/a^,  {-l+a^/ct^),  (2+a^/a^)^ 


Here  a  , ...,a,  are  the  Lagrange  interpolation  coefficients  related 
to  a  second  intersection  between  the  boundary  and  the  mesh  line. 
If  the  mesh  line  in  question  intersects  ^O  in  several  points,  the 
matrix  representing  the  difference  approximation  of  -{'d/'dx.)      along 
this  line  will  be  a  direct  sum  of  several  matrices  of  the  form  (2.3), 
The  matrix 'A.  which  corresponds  to  the  entire  approximation  of 
-{b/bx.  )      is  the  direct  sum  of  the  matrices  introduced  for  the 
individual  mesh  lines  parallel  to  the  vector  e..   Finally,  the 

matrix  A,  which  represents  the  approximation  of  the  entire  problem 

T 
(2.1),  is  the  sum  of  P. A. P.  where  P.  is  a  suitable  permutation 

matrix. 


We  note  that  if  some  irregular  mesh  point  x  is  very  close  to 
the  boundary,  i.e.  some  s  is  quite  close  to  one,  the  ratio  a-, /a 
will  become  very  large.   This  will  give  the  matrix  a  very  large 
diagonal  element  and  the  coefficient  multiplying  g(x*)  in  the 
right-hand  side  will  be  of  the  same  order  of  magnitude.   In 
practice,  we  will  therefore  always  scale  the  rows  of  the  matrix  A, 
making  the  diagonal  elements  equal  to  2n. 


§3.  Stability  of  the  Finite  Difference  Methods 

As  we  saw  in  §2  the  matrix  A  which  corresponds  to  the  full 
difference  approximation  of  problem  (2.1)  has  the  form 


n   ^ 

pIa.p. 

Ill 


5 


where  the  P.  are  suitable  permutation  matrices  and  the  A.  are 
direct  sums  of  matrices  of  the  form  (2.3).   The  original  problem 

(2.1)  has  a  bounded  inverse  in  Lp.   The  analogous  result  is  that 

_2 
the  spectral  norm  of  the  inverse  of  A  is  bounded  by  const  x  h 

To  establish  this  result  we  will  study  the  symmetric  part  of  A. 

In  this  section  we  will  use  the  Euclidean  vector  norm  and  the 

spectral  matrix  norm  exclusively. 

Lemma  1.   Let  the  symmetric  part  of  a  matrix  A  satisfy 
(A  +A^)/2  >  51  ,    5  >  0  . 

Then  A  is  nonsingular  and  |A~  |  <  1/5. 
Proof:   Let  Ax  =  b.   Then 


5x^x  <  x^(A  +  A^)x/2  =  (x^b+b^x)/2  <  |b|-|x|  . 


Thus    |x|    <    I  Ax  1/5   which  proves    the    lemma. 


n 
Lemma  2.      Let  A  =  '^         P^. A.P.    where    the    P.    are   permutation 


matrices.      If 


then 


(A^  +  A^)/2    >   61 


(A  +  A^)/2    >   n5I    . 


The  proof  follows  from  an  elementary  variational  argument.   The 
proof  of  the  next  lemma  is  equally  easy. 

Lemma  3 .   Let  the  matrix  A.  be  the  direct  sum  of  certain 


matrices  B. . .   If 


then 


(B.  .  +  B: . )/2  >  51  ,    for  all  j  , 
J- J    -'-J 


(A^+  A^)/2  >  51 


We  are  now  ready  to  apply  these  lemmas.   Specifically  we  will 
study  matrices  of  the  form  (2.3).   For  technical  reasons  we  will 
assume  that  all  these  matrices  have  an  order  of  at  least  2k-l. 
This  condition  can  again  be  satisfied  for  any  smooth  hO   if  the  mesh 
size  h  is  chosen  small  enough.   We  will  reduce  the  study  of  the 
matrix  (2.3)  to  a  simpler  case  which  corresponds  to  imposing  a 
Neumann  condition  at  one  end  point. 

Lemma  4.   Denote  by  B  the  matrix  defined  by  formula  (2.3), 
Let  Bp  be  a  matrix  of  the  form 


f   1    -1 

-1      2 


V 


-1 


t^/a^,  .  ..,a^/a^,  (-l+ag/a^),  {2+a^/cL^)^ 


10 


and  let  B-,  be  a  matrix  of  the  same  form  generated  by  a  ,.,.,cl. 
Suppose  further  that  the  orders  of  the  matrices  B-, ,  Bp  and  B, 
denoted  by  n-.,    n     and  m  respectively,  satisfy  the  conditions. 


If 


and 


then 


(B^  +  B^)/2  >  51 


(Bp+bJ)/2  >  51 


(B  +B^)/2  >  51 


Proof:   Denote  by  B-,  the  matrix  obtained  from  B-,  by  reversing 
the  order  of  its  rows  and  columns.   The  proof  follows  from  the 
identity, 

x^(B+B^)x/2  =  u^(B^  +  B^)u/2  +  v^(B2+B^)v/2   , 

T  T 

where  u  =  (x-,,...,x   )  and  v   =  (x   ,...,x  ).    This  identity  can 
-L      n-j  n-i      m 

be  verified  straightforwardly.   Hence,  by  our  assumption, 
x^(B+B^)x/2  >  5(u^u  +  v^v)  . 

To  conclude  the  proof  we  only  note  that 

T   ,   T      T 
U  U  +  V  V  >•  XX  . 

T  T 

¥e  will  next  use  the  LDL  factorization  of  S  =  (B^  ^BZ)/2. 

to  verify  that  S  is  positive  definite  and  also  give  a  lower  bound 

for  its  eigenvalues.   We  will  write  S  as  a  block  matrix. 


11 


(3.1) 


'11    ^21 


whe  re  s 


2^=  (0,...,0,  aj^/2a^,  ...,a^/2a^,  -l+a2/2a^)  and 


s  p  =  2+a-,/a  .   Its  block  factorization  takes  the  form 


/  L. 


11 


(J.      Q\  fL^^  i- 


0    d 


where 


■^11 


/  1 


\ 


is  bidiagonalj 


^  =  ^2Al  ^  (0,...,0,a^/2a^,(a^+a^_^)/2a^,... 


(aj^+  ...  +a^)/2a^,  -1+  (a^+  ...  +  a^)/'^o.^) 


and 


d  =  S22  -i^^  =  2  +a-j_/a^  -  {  (a^/2a^)^  +  ... 


((a^+  ...  +  a^  V2a^ 


+  a,)/2aj2+  (-1+  (a^+  .  . .  +  a^  )/2a^  )^] 


By  using  the    fact   that   a    +  . . .  +  a,     =   1,    we   find   that 

(3.2)        d  =   l/a^-[a2  +  (aj^  +  a^_^)2+  ...+  (aj^+  ...  +a2)^}/4a2    . 

Computer  results  show  that  the  rational  function  d  is  strictly 


12 


positive  for  0<_s<_l  and  all  lj<k<6.   For  k  =  7  and  8  it  changes 
sign  in  the  interval.   These  results  can  of  course  also  be  verified 
by  a  tedious  paper  and  pencil  calculation.   We  note  that  d  goes 
to  positive  infinity  when  s  approaches  1  while  the  components  of 
Sp-,  and  S,    remain  bounded.   We  are  now  ready  to  establish  a  lower 
bound  for  the  eigenvalues  of  S. 

Lemma  3.   Let  d  .   denote  the  minimum  of  the  function  d(s) 

defined  by  formula  (3.2).   Then  there  exists  a  strictly  positive 

constant  C,  independent  of  the  mesh  size  h  and  the  region  O,  such 
that 

S  ^  51 
where 

S  =  Cd^.^hV(diajneter  (o))^  . 

Proof:   By  using  the  notations  previously  introduced  in  this 
section,  we  find 

x^Sx  =  x^LDL^x  >  min  (d^^^, 1 ) | L^x | ^  . 

Since  d  =  1  for  s  =  0,  x  Sx  ^  d  |L  x|  .  To  obtain  a  lower  bound 
for  |L  x|  we  will  compute  an  upper  bound  for  | L~  y|.  Partitioning 
the  vector  so  that  y  =  (y  ,y  ),    we  find 


yV^  =  ((y^-y^^)L-J,y^)  . 


Therefore,  if  we  use  the  fact  that  H   has  a  uniformly  bounded  norm, 
we  find 


13 


,-1 


The  norm  of  L-, -,  equals  the  square  root  of  the  reciprocal  of  '^he 


smalle 


st  eigenvalue  of  L-,-,] 


11' 


Now     the   matrix 


T 
^llhl 


/  1 
1 


has  an  order  m  <_   diam  (o)/h.   As  is  easily  checked,  the 

smallest  eigenvalue  of  L-j^]_L-^2_  equals  4  sin  (7r/2  ( 2m+l ))  corresponding 

to  an  eigenvector  with  the  components  cos  (tt  ( j-l/2  )/(2m+l)  ), 

j  =  1, ..o,m.   This  concludes  the  proof. 

By  combining  our  five  lemmas  and  the  results  from  our  compu- 
tation of  d  .  ,  we  obtain,  what  essentially  is  Kreiss'  result. 

Theorem  1.   For  k  _<  6  there  exist  constants  C,  ,  ind'^pendent 
of  h,  such  that, 

(A"^|  <  C^(diajii  o)^  xh~^  . 


14 


§4=  Convergence  and  Asymptotic  Expansions  of  the  Error 

In  this  section,  we  will  prove  the  convergence  of  the  schemes 
introduced  in  §2  and  simultaneously  establish  asymptotic  expansions 
for  the  error.   We  will  concentrate  on  the  case  k  =  6,  which  is 
the  most  accurate  of  the  schemes  known  to  be  stable.   We  will 
assume  throughout  that  the  solution  u(x)  is  sufficiently  smooth. 
We  make  the  Ansatz, 

(4.1)  u^(x)  =  u(x)  +h2e^^^x)  +h^e^^)(x)+  r^(x)  . 

The  functions  e^-^^(x)  and  e^^^(x)  will  be  chosen  as  solutions  of 
Poisson's  equation  in  a  way  which  will  make  the  remainder  r  (x) 
a  term  of  higher  order. 

Asymptotic  expansions  of  this  form  are  basic  for  the  justifi- 
cation of  Richardson  extrapolation  and  deferred  correction  methods. 
They  also  easily  enable  us  to  give  estimates  for  the  rate  at  which 
difference  quotients  of  the  solution  of  the  discrete  problem  u  (x) 

converge   to  the  corresponding  derivatives  of  the  solution  u(x). 

2 
Let  us  denote  by  h  L,  the  difference  operator  which  has  the 

matrix  representation  A,  see  §§2  and  3*   The  linear  system  of 

equations  therefore  has  the  form 

(4.2)  ^^V'^  =  ^"^  ' 

A  component  of  the  right-hand  side  F  ,  which  corresponds  to  a 
regular  mesh  point,  has  the  form  h  f(x)  whereas  a  component,  corre- 
spending  to  an  irregular  mesh  point,  is  a  sum  of  h  f(x)  and  terms 


of  the  form 


g{x^)/a^{s^).      Here  ^^(s^),  0  <  s^  <  1,  is  a  Lagrange 
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polynomial  coefficient  introduced  in  §2.   To  derive  equations  for 
the  error  functions  e^   (x)  and  e^   (x),  we  substitute  the  expres- 
sion (4.1)  into  the  equation  (4.2)  and  expand  the  truncation  error 
in  the  customary  way.   We  first  ignore  the  contributions  from  the 
interpolation  formulas  used  for  the  irregular  mesh  points.   By 
setting  the  fourth  and  sixth  order  terms  of  the  resulting  expres- 
sions equal  to  zero,  we  obtain  the  Poisson  equations 

-Ae^l)  =  (1/12)  ^   (V^x.  )\ 
and 

-Ae^2)  ^  (1/360)  jf-  (V^x^)^u+  (1/12)  jf-  (V^x.)^e^l)  . 

Because  of  the  high  order  accuracy  with  which  the  boundary  condi- 
tion is  approximated,  it  is  appropriate  to  equip  these  equations 
with  zero  Dirichlet  boundary  conditions 

e^^^x)  =  0  ,    e^^)(x)  =  0  ,    x  e   hO    . 

The  functions  e^  '^(x)  and  e^  '^(x)  are,  by  a  standard  result  on 
elliptic  equations,  sufficiently  smooth  functions. 

We  now  derive  a  difference  equation  for  the  remainder  term 
r  (x)o   This  equation  has  the  form 

h\rh  ^  G^  . 

It  is  easy  to  show  that  a  component  of  G  ,  corresponding  to  a 
regular  mesh  point,  is  of  the  form. 
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h^((l/20l60)  J^  (V^x.)^u+ (1/360)  J^  (V^x.)^e^^) 

+  (1/12)  ^  (V^x.  )^e^2b  +  0(h^°)  . 

A  component  of  G  ,  corresponding  to  an  irregular  mesh  point  is  the 
sum  of  a  term  of  this  form  and  of  one  or  more  seventh  order  error 
terms  of  the  Lagrange  interpolation  formulas  (2.2).   The  latter 
terms  are  multiplied  by  factors  l/a  (s.)  which  grow  as 
const/  (1-s.  )  when  s.  -*  1.   It  can,  however,  be  shown,  by  a 
straightforward  calculation,  that  this  increasing  factor  l/a  (s.) 
is  fully  compensated  by  a  decreasing  factor  in  the  error  bound  for 
Lagrange  interpolation,  see  Isaacson  and  Keller  [8  ,  p.  190]. 
These  components  of  G  are  therefore  uniformly  0(h  )  for  all 
sufficiently  smooth  solutions  u(x). 

We  are  now  ready  to  use  Theorem  1  to  obtain  a  bound  for 
r  (x).   It  is  natural  to  work  with  the  norm, 

\\rX   =  (IZ:h"|r^x)|2)'/'  , 


for  which  the  estimate  of  Theorem  1  still  holds.   ¥e  first  estimate 
l2' 


||G  lip.   The  components  of  G   are  0(h  )  for  the  regular  mesh  points 


and  0(h  )  for  the  irregular  mesh  points.   Since  there  are  only  of 
the  order  h"^'^"  ^    irregular  points,  ||G  H^  =  0(h    ).   We  use  Theorem 
1  to  establish. 

Theorem  2.   Let  u  (x)  be  the  solution  of  the  finite  difference 
scheme  with  k  =  6  and  let  u(x)  be  the  sufficiently  smooth  solution 
of  the  differential  equation  (2.1).   Then  there  exist  two 
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sufficiently  smooth  functions  e^  '^(x)  and  e^  '^(x)  such  that 

u^(x)  =  u(x)  +h2e^^^x)  +  h^e^^)(x)+  r^(x)  ,    x  e  O^  , 

where  the  L^  norm  of  r  (x)  Is  0(h-^*''). 

Similar  results  hold  for  smaller  values  of  k.   We  expect  that 
Theorem  2  is  not  sharp.   We  conjecture  that  the  remainder  term 
should  be  of  the  form 

r'^(x)  =  h^e^^^x)  +  0(h^) 

in  the  maximum  norm.   We  are  led  to  this  conjecture  by  results, 
previously  established  by  Bramble  and  Hubbard  [2],    for  the  opera- 
tors of  strictly  positive  type  which  result  when  k  =  1  and  2.   If 
the  estimate  of  Theorem  2  can  be  sharpened  in  this  way,  we  would  be 
justified  in  applying  Richardson  extrapolation  three  times  to 
obtain  a  seventh  order  accurate  method. 
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§5.  Methods  of  Increasing  the  Accuracy 

Richardson  extrapolation  and  deferred  correction  methods  are 
available  to  improve  the  second  order  accuracy  of  the  basic  solu- 
tion u  (x)..   ¥e  will  again  concentrate  on  the  case  k  =  60   We  will 
first  discuss  the  Richardson  extrapolation  method  which  is  simpler 
both  conceptually  and  in  terms  of  its  implementationo 

The  solution  is  first  found  on  a  basic  mesh  O,   and  then  for 

o 
a  sequence  of  refined  meshes  O,  ,  where  h.  =  h  /v., 

1  <  r-,  <  rp  <  o  o  .  o   It  is  very  important  that  the  sequence  (r.] 

grows  slowly  for  multidimensional  problems  since  the  number  of 

variables  grows  rapidly.   The  improved  solution  is  obtained  only  on 

the  intersection  of  the  meshes  O,  .   If  we  require  the  improved 

solution  at  all  points  of  O,   and  use  two  extrapolation  steps,  the 

o 
number  of  mesh  points  on  the  finest  mesh  will  be  at  least  about 

nine  (twenty  seven)  times  larger  in  two  (three)  dimensions.   Core 
storage  can  therefore  easily  be  exhausted  and  less  advantage  can 
also  be  taken  of  the  savings  which  often  can  be  realized  when 
direct  methods  are  used  to  solve  linear  systems  repeatedly. 

If  enough  terms  of  an  asymptotic  error  expansion,  in  even 
powers  of  h,  exist,  we  obtain  improved  solutions  u"?  by  the  recur- 
sion formula, 

-O  '^i 

With  u:  the  restriction  of  u   to  the  intersection  of  the  meshes 

/v  -i  21+2 

O,  .   The  error  uV  -  u  will  be  of  the  order  h  "^   .A  useful 

i 

a  posteriori  error  bound, 

uj-u;^  (uj-uj_^)/(l-  (r.^/r,)^)  , 
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can  also  be  computed,  for  details  see  Bulirsch  and  Stoer  [5]. 

By  using  Theorem  2,    we  can  easily  show  that  two  steps  of 
Richardson  extrapolation  will  give  an  accuracy  of  the  order  Yr'-'^ 
if  we  use  the  scheme  with  k  =  6. 

The  deferred  correction  method  requires  only  one  mesh.   The 
method  has  been  discussed  in  detail  in  a  number  of  papers,  see  for 

example  Pereyra  [13-17] •   Here,  approximations  of  the  leading  terms 

2 
of  the  local  truncation  error  of  the  discrete  operator  h  L,  are 

computed  and  a  corrected  solution  is  then  found  by  solving  an  addi- 
tional system  of  linear  equations  with  the  same  matrix  A  as  before. 
Further  corrections  may  be  obtained  in  s  similar  way. 

We  will  describe  the  variant  of  the  method  which  we  have  used 
in  our  experiments.   In  the  first  step  we  take  into  account  only 
the  first  truncation  error  terms,  resulting  from  the  approximation 
of  (^/^x. )  ,  i  =  1, ...,n,  by  the  three  point  approximations »   We 
know  from  §4  that  these  leading  terms  are 

(5.1)  h^(l/l2)(a/ax.)^u  ,    i  =  l,...,n  o 

We  attempt  to  approximate  them  to  within  0(h  )  by  using  centered 
five  point  differences  of  the  second  order  accurate  solution  u  . 
For  a  periodic  problem  this  procedure  is  very  simple,  but  for  a 
region  with  a  boundary  special  one-sided  differentiation  formulas 
must  be  used  for  the  mesh  points  which  are  within  2h  of  the  bound- 
ary along  a  mesh  line.   One-sided  formulas  can  introduce  additional 
error  terms  for  the  corrected  solution  through  the  special  contri- 
butions to  the  truncation  error  at  the  points  where  these  formulas 
are  employed.   An  additional  correction  step  may  be  justified  by  an 
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asymptotic  error  expansion  of  the  corrected  solution,  but  note  that 
an  unfortunate  choice  of  one-sided  differentation  formulas  would 
lead  to  difficulties  very  similar  to  those  already  discussed  by 
Wasow  [20] . 

This  problem  can  be  avoided  in  a  systematic  way.   Let  x  be 
the  irregular  mesh  point  introduced  in  our  discussion  in  §2.   We 
will  use  high  order  Lagrange  extrapolation,  employing  only  values 
of  the  mesh  function  u  at  x,  x-he.,.,.,  to  obtain  provisional 
values  of  u  at  x+he .  and  x+2he..   The  same  centered  five  point 
differentiation  formula  can  then  be  used  for  all  points  in  O^,  see 
further  discussion  in  §6.   The  expansion  of  the  truncation  error 
which  is  due  to  the  use  of  the  five  point  approximation  of  the 
fourth  derivative  {h/hx. )      will  have  the  same  leading  terms  and 
differ  only  in  a  higher  order  term.   The  order  of  this  higher  order 
term  will  of  course  depend  on  the  degree  of  the  Lagrange  extrapola- 
tion polynomial.   The  Lagrange  polynomial  coefficients  are  the  same 
at  every  point  since  we  use  values  at  mesh  points  only.   The 
approximation  of  the  expressions  (5.1)  are  added  to  the  original 
data  F   and  the  linear  system  of  equations  is  solved  a  second  time. 

In  a  second  correction  step 

(5«2)     h^(l/l2)(VSx.)^u+h^(l/360)(V^x.)^u  ,    i  =  1,  ...,n  , 


is  approximated  by  a  centered  seven  point  formula  with  an  error 
which  is  0(h  )  for  a  sufficiently  smooth  function.   We  thus  use 
the  once  corrected  solution  and  high  order  extrapolated  values 
thereof  in  a  way  very  similar  to  the  previous  step  to  obtain  a  new 
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right-hand  side  and  a  second  corrected  solution,  see  further 
discussion  in  §6. 

Our  error  bounds  for  the  deferred  correction  method  are 
rather  weak.   When  we  estimate  the  truncation  error  due  to  the  dis- 
cretization of  the  expression  (5.I),  we  find  that  the  three  first 
terms  of  the  expansion  given  in  Theorem  2  give  a  contribution  of 

the  order  h  .   Since  the  operator  h  L,  has  an  inverse  bounded  by 

-2  4 

const  h   they  contribute  a  term  of  the  order  h  to  the  error  of 

the  corrected  solution.   In  contrast  the  undivided  differences  of 
the  remainder  term  of  r   create  difficulties.  Since  undivided 
difference  operators  are  bounded,  independently  of  h,  the  contribu- 
tions of  r   to  the  truncation  error  and  the  error  of  the  corrected 
solution  are  bounded  by  h^  -^  and  h  -^ ,    respectively.   In  order  to 
prove  a  result  as  strong  as  that  for  the  Richardson  extrapolation 
method  this  loss  of  two  powers  of  h  must  be  eliminated.   This  would 
be  achieved  if  we  were  able  to  give  as  sharp  a  bound  for  the  norm 
of  the  second  order  divided  differences  of  the  solution  as  for  the 
norm  of  the  solution  itself.   The  analogue  of  this  desired  estimate 
holds  for  second  order  elliptic  equations  on  regions  with  suffi- 
ciently smooth  boundaries.   We  have  not  been  able  to  obtain  this 
result  in  the  discrete  case.   A  modification  of  the  argument  of  §3 
leads  to  an  improved  bound  for  divided  differences  of  the  first 
order.   This  proves  that  at  most  one  power  of  h  can  be  lost  in  each 
connection  step.   For  numerical  evidence  see  §7. 
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§6.  The  Capacitance  Matrix  Method 

All  our  experiments  have  been  carried  out  for  regions  in  the 
plane  and  we  will  therefore  confine  our  discussion  to  that  case. 
We  have  used  a  modification  of  the  capacitance,  or  imbedding, 
method  which  was  developed  by  Proskurowski  and  Widlund  [18]  to 
solve  our  linear  systems  of  equations.   We  refer  to  that  paper  for 
a  detailed  discussion  of  the  method.   Here  we  will  confine  our- 
selves to  a  few  brief  remarks  on  the  method  concentrating  on  the 
changes  required  by  the  deferred  correction  algorithm. 

A  main  part  of  any  capacitance  matrix  program  is  a  fast 
Poisson  solver  on  a  region  for  which  separation  of  the  variables 
can  be  applied.   Our  subroutine,  SOLVE,  implements  a  Fourier- 
Toeplitz  method  on  an  infinite  parallel  strip  with  periodic 
boundary  conditions  in  one  direction,  see  Proskurowski  and  Wldlimd 
[18,  Section  6]  and  Fischer,  Golub,  Hald,  Leiva  and  Widlund  [6  ], 
Our  region  O  is  imbedded  in  a  rectangular  subset  of  this  strip.  The 
fast  solver  requires  of  the  order  mn  log^n  operations  for  the  exact 
solution  of  the  five  point  discrete  Poisson  equation.   Here  n,  the 
number  of  mesh  points  across  the  strip,  is  preferably  a  power  of  two 
and  m  is  the  number  of  mesh  points  used  along  the  strip.   We  will 
see  below  that  it  is  convenient  to  place  the  region  O  inside  a 
centered  subset,  of  size  (m-6)x  (n-6),  of  the  set  of  mxn  mesh 
points  which  is  used  by  SOLVE. 

An  extended  system  of  linear  equations  with  a  matrix 

r\j  rn 

A  =  B+  UZ 
is  solved.  The  matrix  B  corresponds  to  the  five  point  formula  on  the 
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strip  while  A  contains  our  matrix  A,  see  §§2  and  3,  as  a  principal 
minor.   The  matrices  U  and  Z  are  sparse  and  have  p  columns  where  p 
is  the  number  of  irregular  mesh  points.   The  matrix  U  is  chosen  so 
that  UVj  V  any  p-vector,  is  an  extension  by  zero  of  the  corre- 
sponding mesh  function  v  defined  only  on  the  set  of  irregular  mesh 

T 
points.   The  matrix  Z  is  thus  a  compact  representation  of  the 

matrix  A-B  from  which  zero  rows  have  been  eliminated.   A  change  of 


the  approximation  of  the  boundary  condition  involves  a  change  of 
the  matrix  Z.   The  right-hand  side  F   of  our  original  system  of 
equations  is  extended,  in  an  arbitrary  way,  to  the  complement  of  Q,  . 
The  matrix  A  is  constructed  in  such  a  way  that  the  restriction  of 
the  solution  of  the  extended  system  to  the  set  O,,^  is  the  solution 
of  our  original  system  of  equations. 

There  are  two  main  parts  of  our  capacitance  matrix  program. 
We  execute  the  first  one  only  once  for  a  particular  choice  of  h 
(a  mesh  size),  k  (a  member  of  our  family  of  difference  schemes)  and 
a  region  Oo   In  this  first  part  a  px  p  nonsymmetric  dense  capaci- 
tance matrix  C  is  computed  at  an  expense  of  one  call  of  the  sub- 

2 
routine  SOLVE  and  of  the  order  p   additional  operations.   A  solu- 
tion for  a  specific  set  of  data,  which  is  accomplished  in  the 
second  part,  requires  essentially  only  two  calls  of  the  subroutine 
SOLVE  and  the  solution  of  a  capacitance  matrix  system  of  equations 
of  the  form  C\x   =   g.   In  our  implementation  the  capacitance  matrix 
C  is  very  well  conditioned  and  this  equation  can  therefore  be 

solved  accurately  by  a  conjugate  gradient  method  at  an  expense 

2 
of  the  order  p   operations.   We  have  however  chosen  to  use 

Gaussian  elimination.   The  matrix  C  is  factored  only  once,  at  an 
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expense  of  the  order  p   operations,  and  the  factors  are  then  stored 

and  used  for  any  additional  set  of  data.   Any  subsequent  problem 

o 
therefore  requires  only  of  the  order  (mn  logpn  + p  )  operations. 

The  method  is  numerically  very  stable  and  the  linear  system  of 
equations  is  solved  very  accurately. 

Two  rectangular  arrays  of  the  dimension  mxn  are  used  to 
store  the  data  and  the  solution.   The  first  array  initially  con- 
tains the  original  data  F  ,  arbitrarily  extended  to  the  complement 
of  O,  .   The  second  order  accurate  solution  u  is  computed  and 
stored  in  the  second  array.   This  solution  is  then  extended  to 


§5.   A  first  contribution  to  the  modified  right-hand  side  of  the 
equation  is  then  computed  by  using  a  five  point  numerical  differen- 
tiation formula  on  all  mesh  lines  parallel  with  the  x-,-axis.   The 
resulting  mesh  function  is  added  to  F  ,  the  content  of  the  first 
array.   This  process  is  now  repeated  in  the  other  direction.   ¥e 
thus  extrapolate  u  (x)  in  the  x^-direction  to  the  appropriate  ex- 


tion  to  obtain  the  final  contribution  to  the  new  right-hand  side. 
We  note  that  we  can  simplify  the  programming  by  using  the  numerical 
differentiation  formula  over  the  entire  rectangular  region  since 
the  restriction  of  the  solution  on  the  strip  to  the  set  O,  is 
independent  of  the  values  of  the  data  outside  n^.^.   The  second  part 
of  the  capacitance  matrix  solver  is  now  used,  with  the  new  right- 
hand  side,  to  produce  a  fourth  order  accurate  solution.   It  is 
stored  in  the  second  array  which  also  serves  as  work  space  during 
this  part  of  the  calculation.   The  final  corrected  solution  is 
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computed  similarly.   The  original  data  F   is  read  into  the  first 
array  and  approximations  to  the  expressions  in  formula  (5.2)  are 
added.   In  this  step  seven  point  differentiation  formulas  are  used. 
We  note  that  since  we  placed  O,  inside  a  rectangle,  leaving  three 
extra  mesh  lines  on  all  sides,  we  can  carry  out  all  the  necessary 
extrapolations  while  using  only  the  storage  locations  provided  for 
in  the  second  mxn  array.   This  admittedly  introduces  an  additional 
constraint  on  the  choice  of  mesh  size  for  certain  nonconvex  regions 
but  this  aspect  of  the  implementation  of  our  method  can  of  course 
easily  be  changed.   The  extrapolation  and  numerical  differentiation 
steps  are  very  straightforward  and  require  very  little  computer 
time,  see  §7. 
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§7 •  Numerical  Experiments 

A  FORTRAN  program  incorporating  the  ideas  of  this  paper  was 
prepared  and  run  in  single  precision  (between  l4  and  I5  decimal 
digits)  on  a  CDC  760O  computer  at  the  Lawrence  Berkeley  Laboratory- 
using  a  RUN  76  compiler.   We  report  on  experiments  using  second 
and  sixth  order  Lagrange  interpolation  formulas,  k  =  2  and  6,  for 
the  irregular  mesh  points,  see  §2o   In  all  our  experiments  the 
region  was  a  circle  of  radius  one  centered  at  the  origin  and  the 
mesh  size  was  h  =  I/23 .   There  were  I653  mesh  points  of  which  128 
were  irregular  and  the  region  was  imbedded  in  a  6k  xGk   mesh. 

By  e   and  £„  we  denote  the  maximum  and  Lp  norms  of  the 


error,  i.e. 


and 


max  |u  (x)  -  u(x) |  , 

{(1/N)  y—    |u^(x)-u(x)|2]^/^ 


In  Table  1,  we  report  on  the  solution  of 

-Au(x)  =  2  sin  (x-j^+x^) 

with  boundary  values  and  exact  solution  equal  to  u(x)  =sin  (x^+  x  ) 
This  is  a  problem  with  a  very  smooth  solution  and  served  basically 
as  a  test  that  the  program  and  algorithm  really  worked.   We  note 
that  we  obtain  close  to  full  word  accuracy. 
The  next  problem,  see  Table  2,  was 

-Au(x)  =  53  sin  (2x-^  -  7x2) 
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with  the  boundary  values  and  exact  solution  equal  to 

u(x)  =  sin  (2X-,  -  7Xp ) .   This  problem  is  more  difficult  than  the 

first  since  the  solution  is  more  oscillatory.   We  tried  sixth  and 

second  order  interpolation  at  the  irregular  mesh  point.   According 

to  results  of  Bramble  and  Hubbard  [2  ]  there  is  an  expansion  of  the 

form 

u'^(x)  =  u(x)  +  h^e^^^x)  +  0(h^) 

when  second  order  interpolation,  k  =  2,  is  used.   We  note  that  the 
first  correction  step  gives  a  smaller  improvement  in  the  case  k  =  2 
than  when  k  =  6  and  that  the  second  correction  step  gives  no 
improvement  for  k  =  2.   This  experiment  thus  confirms  the  observa- 
tions of  Wasow  [20],  Pereyra  [13]  and  others  on  the  importance  of 
the  existence  of  asjrmptotic  error  expansions.   We  also  note  that 
the  two  second  order  methods,  obtained  before  the  correction  steps, 
perform  equally  well. 

A  final  series  of  experiments  were  carried  out  to  study  the 
effects  of  lack  of  smoothness  of  the  solutions.  The  problems  had 
the  form 

f-2£{jl-l){x^+x^)^'^,       if      x^+X2>0 
-Au  =  I 

I  0  ,       otherwise 

with  the  boundary  values  and  exact  solution  equal  to 

r(x^+  x^)^  ,   if   (x^+  x^)  >  0 
u(x)  = 

[_    0     ,   otherwise  . 

We  tried  i  =  2,  4  and  6.   The  solution  then  has  a  jump  discontinuity 
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in  derivatives  of  order  £.   The  results  are  given  in  Table  3.   The 
performance  of  the  method  with  k  =  2,  H,   =  6,    is  consistent  with 
our  previous  observations.   For  k  =  6  and  with  i  =  2,  4  it  appears 
as  if  a  ^-th  order  accurate  method  is  obtained  for  these  solutions 
which  have  a  jump  in  the  i-th  derivatives.   Care  must  of  course  be 
exercised  when  trying  to  draw  such  conclusions  from  our  very  limi- 
ted experimental  evidence.   We  feel  however  that  our  results  are 
encouraging.   We  note  that  when  the  solutions  fail  to  be  smooth 
enough  the  corrections  do  not  destroy  the  accuracy  obtained  in  the 
previous  steps. 

The  total  CPU-time  for  a  problem  with  k  =  6  was  10.28 
seconds.   The  first  part  of  the  capacitance  matrix  program,  see  §6j 
computed  the  second  order  accurate  solution  u  (x)  in  8.77  seconds. 
The  first  correction  required  an  additional  0.66  seconds  and  the 
second  correction  took  an  additional  O.85  seconds.   In  the  correc- 
tion steps  the  extrapolation  to  exterior  mesh  points  and  the 
differentiation  steps  required  less  than  10^  of  the  time.   The 
execution  time  could  be  reduced  by  optimizing  our  program  and  by 
changing  to  a  faster  compiler. 


Correction 

0 

1 

2 

e      ,    k  =  6 

CO  ' 

1.9  xlO"5 

1.0  XlO"^ 

5.6x10-1^ 

62   ,    k  =  6 

1.0  xlO"5 

5.4x10-10 

3.4  X  10-12 

Table  1 

Lp-  and  maximum-norm  errors  for  a  problem  with  the 
solution  u(x)  =  sin  (x^+X2).   Sixth  order  inter- 
polation used  at  the  boundary  points. 
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Correction 

0 

1 

2 

-00'    ^   =   2 

8.8  XlO"-^ 

1.3  xlO"^ 

1.4  xio"^ 

eg   ,    k  =  2 

4.7  xio"^ 

3.3  xio"^ 

3.4  xio"^ 

^oo'    ^  =  6 

9.2  xio"^ 

5.3  xlo"5 

1.3  xio"^ 

£2    .    k  -  6 

4.8  xio"^ 

2.8  X10"5 

3.4  xio"^ 

Table  2 

Lp-  and  maximum-norm  errors  for  a  problem  with  the 
solution  u(x)  =  sin  (2x-,  -  7Xp).  Second  and  sixth 
order  interpolation  are  used. 


Correction 

0 

1 

2 

B^,     ^    =   2,    k   = 

6 

9.9x10"^ 

9.9  xio"^ 

9.9  xio"^ 

B^,      ^     =4,     k    = 

6 

1.2  xlO"^ 

4.8  xio"^ 

4„8  xio"^ 

B^,     i    =   6,    k   = 

6 

7.4  xlo~^ 

9.4  xio""^ 

7«7  xio"^ 

B^,      i     =    6,     k    = 

2 

6.7  xio"^ 

1.5  xio"^ 

1.5  xio"^ 

Table  3 

Maximum-norm  error  for  a  problem  with  the  solution 
u(x)  -  (x-^  +  Xp)^,  if  (x-^+Xg)  >  0,  u(x)  =  0  other- 
wise.  Sixth  and  second  order  interpolation  are  used, 
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